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Abstract 


The most important purpose in location problems is usually to locate some 
facilities and allocate the demands of nodes so that the total transportation 
cost of the network is minimized. However, in real networks, there are some 
other influencing factors, aside from the transportation costs, for determin- 
ing the allocation mode. In this paper, a minimum information approach is 
applied to the capacitated p-median problem to estimate the most likely allo- 
cation solution based on some prior probabilities. Indeed, the most probable 
solution is achieved through minimizing a log-based objective function, while 
the total transportation cost should be less than or equal to a predetermined 
budget. The problem is solved by using a decomposition method combined 
with the Karush-Kuhn—Tucker optimality conditions, and some numerical 
examples are provided to verify the added value of the proposed model and 
solution approach 
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1 Introduction 


The p-median problem is one of the known subjects in the field of locating 
facilities with a minisum objective function, where its purpose is to locate 
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p facilities and allocate the demands of nodes so that the total demand- 
weighted distance is minimized. ReVelle and Swain [26] introduced the first 
integer linear programming formulation for the p-median, which has played an 
important role in solving the problem by applying such exact methods as the 
branch and bound [24] and Lagrangian dual relaxation [5]. Hakimi [18] proved 
that the p-median problem on the general networks is NP-hard; however, it 
has at least one nodal solution. Kariv and Hakimi [20] suggested an O(p?n?) 
time algorithm to solve the p-median problem on the tree networks. Then, 
Tamir [30] improved the time complexity on the tree networks to O(pn?). 
For a complete bibliography on the p-median problem; see [25]. 

In the classical p-median model, which is sometimes referred to as the 
uncapacitated p-median (UC'pM) problem, there is no constraint on the ca- 
pacities of facilities; see [17]. In the case of limited capacities, the problem is 
usually called the capacitated p-median (CpM) problem. Both the UCpM 
and C'pM problems are NP-hard, although considering capacity limitations 
makes the latter problem more difficult; see [14]. A column generation ap- 
proach, along with a Lagrangian/surrogate relaxation technique has been 
applied by Lorena and Senne [21] for solving the CpM problem. Ceselli and 
Righini [8] presented a branch and price algorithm that utilizes a column 
generation, heuristics, and branch-and-bound to determine the optimal so- 
lutions. Such Meta-heuristic methods as a Hybrid scatter search and path 
relinking algorithm [10], variable neighborhood search [12], neural network 
[27], and genetic algorithm [23] have been recently applied to solve the CpM 
problem. 

When there is a lack of information in the network, the maximum entropy 
(ME) principle or its extended concept, the minimum information theory, 
is used to consider the solution with the most unbiased probability distribu- 
tion; see [22]. In the ME models, there is no assumption about the missing 
information, and all possible states consistent with the system conditions are 
considered equally; then the most likely one is selected; see [9, 19]. The ap- 
plication of ME to the location problems was employed by Teye, Bell, and 
Bliemer [32, 33] to solve a multiuser inter-modal terminal (IMT) location 
problem wherein there is no priority to use an IMT. Based on the ME 
principle, the most unbiased probability distribution was achieved through 
considering all possible states of [MT usage and selecting the most likely 
one consistent with the problem constraints. Teye, Bell, and Bliemer [31] 
also proposed an MF approach for locating competitive multiuser freight 
facilities in general, and inland multiuser [MTs, in particular, when the 
multiple users have the choice whether to use a facility or not. 

The ME approach to the origin-destination (O — D) matrix estimation 
problems was studied by Van Zuylen and Willumsen [34] to determine the 
most probable O — D demands consistent with the available traffic counts. 
Abareshi, Zaferanieh, and Keramati [2] introduced an ME path flow estima- 
tor for disaggregated flows between O — D pairs with a prespecified level for 
each disaggregation. 
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When, in the network, there is some prior information, but insufficient, 
it is reasonable to choose a solution that adds the least extra information to 
the available knowledge, which results in the most possible unbiased prob- 
ability distribution. The minimum information (WJ) approach provides an 
extended measure of the likelihood of a certain macro-state on the existence 
of some appropriate micro-state space; see [28]. Abareshi and Zaferanieh [1] 
proposed a new bi-level p-median problem in which the total cost of locating 
facilities and serving demands was minimized through the upper level while 
the MI approach was applied in the lower level to determine the most unbi- 
ased allocation solution based on some prior partial information. The prior 
information was given in terms of the probabilities of serving the demands of 
client nodes by different facilities, which might be determined by considering 
some attributes such as distance, locality, and geographical features, taken 
into account by multiple attributes decision making procedures; see [1, 29]. 

In this paper, considering some prior probabilities for serving the demands 
of nodes, we attempt to locate p capacitated facilities and determine the most 
probable allocation solution, while the total transportation cost should be 
less than or equal to a predetermined budget. The model is formulated as a 
mixed-integer nonlinear programming (MINLIP) problem, which is known 
to be NP-hard. Applying such alternative methods as the Lagrangian dual 
approach or Benders decomposition algorithm to separate the model into 
some easier subproblems, the solution could be iteratively obtained. 

The generalized Benders decomposition (GBD) algorithm has been sug- 
gested by Geoffrion [16] as a solution procedure for certain nonlinear pro- 
gramming (NLP) and MINLP problems. The decomposition procedure 
was first developed by Benders [6] for the solution of mixed-variable pro- 
gramming problems. However, some restrictions regarding the convexity and 
other properties of the involved functions were considered. Floudas, Ag- 
garwal, and Ciric [13] studied this technique and revealed its potential for 
applying in a chemical process design. They also proposed a computational 
implementation to reach the global optimal solution for nonconvex NLP and 
MINCTP problems. This claim was sustained by solving different examples, 
though it has not been proven mathematically. 

Bagajewicz and Manousiouthakis [3] investigated the importance of two 
properties defined by Geoffrion [16], called property (P) and L-dual-adequacy, 
for the proper application of the GBD algorithm. They solved several ex- 
amples satisfying or violating these properties and analyzed the obtained 
solutions. Through examples and mathematical analysis, they showed that 
dual gaps may prevent this procedure from converging to the global solution, 
and therefore only bounds on the global optimum might be obtained. Here, 
we apply the GBD algorithm studied in [3] to solve the proposed MINLP 
model of the p-median problem in the MI approach. 

The remainder of this paper is organized as follows: In Section 2, the 
model of the CpM problem in the MI approach is introduced. Using the 
scheme proposed by Bagajewicz and Manousiouthakis [3], the GBD algo- 
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Table 1: The notation used in the model 


Notation Definition 


N {v1,...,Un}, the set of nodes 
E The set of links 
I The index set of client nodes 
J The index set of candidate points for establishing facilities 
Yj The binary variable representing whether node v;,7 € J is used to locate a facility or not 
Cj The capacity of the facility at node v;,7 € J 
Wi The total demand of the client node v;, i € I 
Dig The probability of serving the demand of node v;,7 € I by the facility at node v;,j € J 
Lij The decision variable representing the amount of demand of node v;,i € I 
provided by the facility at node v;,7 € J 
di; Per unit cost to serve the demand of node v;,i € I by the facility at node vj, 7 € J 
The total available budget 


rithm and the Karush—Kuhn—Tucker (kK KT) optimality conditions are ap- 
plied in Sections 3 and 4 whereby the problem is reduced to easier subprob- 
lems. The feasibility and optimality constraints are added iteratively and the 
upper and lower bounds are updated. The convergence of the GBD algorithm 
on the proposed problem is also established in Section 4. To investigate the 
added value of the proposed model, some numerical examples are provided 
in Section 5, while the efficiency of the implemented algorithm is compared 
with some existing MI NIP solvers by applying the performance profile test 
proposed by Dolan and More [11]. The summary and conclusions are given 
in the last section. 


2 CpM problem with the most likely allocation solution 


In this section, the CpM model applying the MI approach, to determine 
the most probable allocation solution, is proposed. Consider a network G = 
(N, £), where the other frequently used notations are listed in Table 1. 

To apply the concept of MI to the CpM problem, consider several de- 
mand nodes v;, 7 € I with w; clients that should be served by the facilities 
at candidate nodes vj, j € J. Also, let xj; be the amount of clients of the 
demand node v; that are served by the facility at node vj, where its corre- 
sponding probability is given by p,;. Then, using the multinomial distribu- 
tion, the probability of observing an allocation solution with the components 
Xi; is given as follows (see [1]): 


pe) = WG Ll) 


Following the scheme presented by Abareshi and Zaferanieh [1], by taking 
the logarithm, using Stirling’s approximation, removing constant values, and 
reducing the problem to the minimum case, the model of the minimum infor- 
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mation capacitated p-median (MICpM) problem is introduced as follows: 


MICpM ante Zo(“,y) = » >) xij (Ina,;; — 1—Inp,;), (1) 
jEeJ iel 

> =P, (2) 
jed 
Stig = 0; for all i € J, (3) 
jed 
¥ 23 < ye; for all j € J, (4) 
i€l 
> eee < B, (5) 
jeJ ier 
yj; € {0,1}, x7 > 0 for alli € I,j € J. (6) 


Constraint (2) assures that the number of established facilities is equal 
to p and constraints (3) practically guarantee that the demands of all nodes 
vu; are supplied. Constraints (4) make all open facilities serve the demands 
less than or equal to their capacities. In addition, constraint (5) states that 
the total spent cost does not exceed the available budget B. Due to the ex- 
isting of nonlinear terms and mixed-integer variables, the MIC'pM problem 
(1)-(6) is NP-hard; see [14]. In the following, the GBD algorithm is intro- 
duced whereby the problem is reduced to smaller subproblems to which the 
constraints are added iteratively to update the upper and lower bounds. 


3 Generalized Benders decomposition algorithm 


Here, the GBD algorithm studied by Bagajewicz and Manousiouthakis [3], 
for general MI NLP is stated, while its application to the MIC'pM problem 
(1)—(6) is examined through the next section. Consider the following problem 


min F(x, y), (7) 
Ly 

s.t. G(a,y) <9, 

zExX, yeEY. 


The projection of the above problem on Y was proposed by Geoffrion [16] 
as 


min v(y), (8) 


s.t. u(y) = min{ F(x, y) : st. G(x,y) < O}, 
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where V = {y: G(a,y) < 0 for some « € X}. Introducing the following 
two problems, Geoffrion [16] proposed a decomposition of (8) when X is a 
convex set, and the functions F(a, y) and G(a,y) are convex with respect to 
the variable «x: 


[1.] Primal problem: 


where y is an arbitrary but fixed point in Y. 
[2.] Master problem: 
P : t 
min{max{min F(x, y) + w' G(s, y)}}, (10) 
s.t. min{'G(x, y)} <0 forallA€ A, 


where A= {\E€R™:A>0, Soy, Ai = 1} (m is the size of vector G). 


Bagajewicz and Manousiouthakis [3] pointed out that the master problem 
is equivalent to the projection (8), in the case, where X is a convex set 
and F'(a,y) and G(a, y) are convex functions with respect to the variable z. 
In addition, they showed that the master problem is also equivalent to the 
following problem: 


min yo, 11) 
YsYo 
s.t. L*(y,u) = min{ F(a, y) +u'G(z,y)}<yo for allu>0, 12) 
“rE 
L..(y,) = min{A'G(a, y)} <0 forallAE A. 13) 
Geoffrion [16] suggested to solve a relaxed version of problem (11)—(13) 


in which all but a few constraints are ignored. Indeed, in each iteration, 
the constraints corresponding to a subset of vectors u and A, instead of all 
of them, are considered. Due to ignoring some of constraints, the obtained 
problem is called the relaxed master problem; see [3]. Since constraints are 
continuously added, the optimal values of this problem form a monotone 
nondecreasing sequence as the lower bounds to problem (7). Next, the steps 
of the GBD algorithm in the general form are represented. 
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3.1 General form of the GBD algorithm 


Step 1. Let a point 7 € YN V be available. Solve the primal problem (9) 
and obtain the optimal solution «* and the optimal multiplier vector 
u*. Put the counters Kf = 1, K' = 0. Set UB = F(x*,y). Select 
a tolerance € > 0 and set uw“) = u*. Finally, determine the function 
E*(y,ul). 


Step 2. Solve globally the current relaxed master problem: 
min yo, 
Y¥,Yo 
s.t. L*(y, ul*s)) <yo, ki =1,...,K’, (optimality constraint) 
L(y, A%)) <0, ko =1,...,K*, (feasibility constraint). 


Let (§, 90) be the globally optimal solution. Set LB = jo. If UB < 
DB +e, then terminate. 


Step 3. Solve globally the primal problem (9) using % = 9. 


Step 3.a. Feasible primal problem: Solve the primal problem (9) by using 
gy. If vg) < LB +e, then terminate. Otherwise, determine the 
optimal multiplier vector u*, set K/ = Kf +1 and uk?) = ys. If 
u(y) < UB, then set UB = v(y). Finally, determine the function 
L*(y,u)) and return to Step 2. 

Step 3.b. Infeasible primal problem: Determine a set of values of A* € A 
satisfying minzex{(A*)*G(z, y)} > 0. Set K’ = K* +1 and 
MK) = \*. Determine the function L,(y, A“). Return to Step 
2 


Remark 1. To determine the value of \* in Step 3.6, the following problem 
should be solved (see [3]): 


mina, 
Lia 
s.t. G(x,9) -—al <0, 
rExX, a€R, 


where 1 = (1,1,...,1)*. 


Remark 2. If the functions F(x, y) and G(x, y) are separable with respect 
to « and y and convex with respect to x, the global optimal solutions for x 
within constraints (12) and (13) are obtained independently of y; see [3]. In 
this case, functions L*(y, ulk’)) and L,(y, A“')) are defined in explicit forms. 

Geoffrion stated some conditions ([16, Theorem 2.4] for the finite discrete 
set Y and [16, Theorem 2.5] for the infinite cardinality set Y) for problem (7), 
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where the lowest solution of the primal and the global solution of the master 
would approach one another, and consequently, the global optimum of the 
overall problem would be provided within a prespecified tolerance. However, 
when the convexity in x does not hold, dual gaps may exist. Nevertheless, the 
global solution to the master problem will still provide a valid lower bound to 
problem (7). Next, an application of the GBD algorithm on problem (1)-—(6) 
is investigated. 


4 Implementing the GBD algorithm on the MICpM 
problem 


To apply the proposed GBD algorithm, we reformulate problem (1)-(6) as 
follows: 


min Zo( x,y) = Sau Inz;; —1—Inp,;), 


jEJ iel 
So aig < yjcj for all j € J, (14) 
icl 
ctEX={aj>0: > ay =u; foralliel, S> >  dijaiy < Bh, 
jed jeJd iel 
yeY=t{y: y,; is binary and Soy; = 


jEed 


Note that the set X as well as functions Zo(x,y) and G;(x,y) = 
ier Tig — ¥jCj, for all j € J, are convex with respect to x. Following 
the steps of Algorithm 3.1, the implementation of the GBD algorithm to 
solve the MICpM problem is summarized as below. 


4.1 GBD algorithm on the MICpM problem 


Step 1 Select a tolerance € > 0 and a feasible vector y € Y for which there 
exists a feasible solution x € X satisfying in (14). 


Step l.a. Solve the corresponding primal problem (15)—(16) 


min Zo(2, ¥) => >> 255( lng; —1—Inp;,), (15) 
jeJ ie] 
Ya Se for all j € J, (16) 


tel 
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and obtain the optimal solution «* and the optimal multipliers 
u;, j € J corresponding to constraints (16). Set Kf =1,K'=0, 
UB = Zo(a*,¥) and vector uw) = u*. 


Step 1.b. Determine the optimal solution < corresponding to the following 
subproblem and set 2") = #: 


min >) >) tig (In ay — 1 —Inpig) + Sou) wg. (17) 


Jed i€l jEed tel 


Step 2 Solve globally the current relaxed master problem: 


18 
min, vo (18) 
st. L*(y,u(*) aya laces =1—1ing;7) (19) 
jEeJ iel 
+ Su) (57 oft) —¢jyj)< yo, hm =1,...,K!, 
JET tel 
La (y, A) = > NMS a — c7y;) <0, ka =1,...,K4 AEA 
jEed wel 


Note that in each iteration, the variables wu *1),2(*), and a/(*2) are 
specified fixed values that have been obtained during the previous steps 
of the algorithm; therefore, problem (18) is a linear mixed-binary model, 
which could be globally solved by such existing methods as branch and 
bound. Let (%,%o) be the global optimal solution. Set LB = Gj. If 
UB<LIB+e, then terminate. 


Step 3. Solve globally the following problem to examine the feasibility of the 
primal problem (15)—(16) corresponding to 7 = @. 


mina, 
So aig — 9; a <0 forall je J, (20) 
rExX, aER. 


Step 3.a. If a < 0, then the primal problem (15)—(16) is feasible. Solve the 
primal problem for the given vector ¥ to find the optimal solution 
x*. If Zo(a*,y) < LB+e, then terminate. Otherwise, determine 
the optimal multiplier vector u* corresponding to constraints (16), 
and set Kf = Kf +1 and ul") = u*. If Zo(a*,y) < UB, then 
put UB = Z)(a*,y). Finally, return to Step 1.b. 


Step 3.b. If a > 0, then the primal problem (15)—(16) is infeasible. Deter- 
mine the optimal multipliers \j corresponding to constraints (20). 
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Set Ki = Ki +1 and \(") = )*. To find the optimal solution 4, 
solve the following model: 


min SA) 


jes ie€l 
Put 2K’) = & and go to Step 2. 


The purpose of the algorithm is to minimize the gap between the upper and 
lower bounds. Next, the convergence conditions stated by Geoffrion [16], are 
verified for problem (1)—(6). 


Theorem 1. The GBD Algorithm 4.1 for the MICpM problem (1)-(6) 
terminates in a finite number of steps, for any given € > 0. 


Proof. Since Y = {y : y; is binary and >7,.7 yj = p} is a finite discrete set, 
the convergence of Algorithm 4.1 should be established via [16, Theorem 2.4]. 
Indeed, the following conditions must hold: 


i. X is a nonempty convex set and functions Gj(x,y) = doje, Lig — YC; 
are convex on X for each fixed y € Y. Further, the set Z, = {z = 
(21,---,2 |) : Gj(x,y) < 2; for some x € X} is closed for each fixed 
yeY. 


ii. Functions Zo(z, y) and G(x, y) are convex on X for each fixed y € Y. 
In addition, for each fixed y¥ € YNV, either the optimal value Zo(zx, 7) is 
infinite, or the primal problem (15)—(16) possesses an optimal multiplier 
vector. 


Part i. Obviously, the set X = {x;; > 0: Pes Lig = wi, Vi € 
I, Mies Vier Uijtij < B} is nonempty (xj = wi, for allie IJ and xj; = 
0 for 2 # 7 is a feasible solution), closed and bounded. Noting the continuity 
of functions Gj (x,y) on X, the closeness of Z, is implied; see [16]. 

Part ii. The convexity of functions Zo(2,y) = dijey Mier Vij (IM 2iy — 1 — 
Inpij) and G;(z,y) = jie, Lig — yjcj with respect to x for each fixed y 
is straightforward by the definition. Following the saddle point optimality 
conditions stated in [4], the existence of optimal multiplier vector is inferred; 
see [15]. 


The overall scheme of the algorithm is represented in Figure 1. 


4.2 Solving the primal problem 


The primal problem (15)—(16) is a convex nonlinear problem, which can be 
solved by applying the K KT optimality conditions. Let y be given. Con- 
sider the Lagrangian dual vectors u;,7;, and 7 corresponding to constraints 
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>| Solve primal problem (4.2) 


| Update the upper bound 


Solve Problem (4.4) 


Solve master problem (4.5) " No Yes 


Update the lower bound 


| Return the solution 


Solve problem (4.6) 


Yes 


(<0?) No Solve Problem (4.7) 


Figure 1: The overall framework of the proposed algorithm 


(22), (23), and (24). Then Lemma 1 expresses the relationship between the 
primal and dual variables: 


min Zo(z, y) = S- » x43 (In i; — 1 — Inpi;), (21) 


jeJ i€l 
Say = ChY for all ged, (22) 
wel 
Stig = 0; for alli € I, (23) 
JET 
DD dijtiy $B, (24) 
jeJ i€l 
aij >0 forallie l,j € J. (25) 


Lemma 1. The optimal solution of problem (21)—(25) is obtained by 
tig = pie UItvnds), for all j € J,i € I, where u;,7, and 7 are the 
optimal multiplier vectors corresponding to constraints (22), (23), and (24), 
respectively. 


Proof. The K KT optimality conditions of model (21)—(25) with respect to 
x obtain in the following expressions: 


In Lig = In piz + Uj + Vi + ndij > 0, (26) 
xij (Inv, Inpyy +uj+y%4 ndi;) = 0. (27) 
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The above conditions state that one of the following cases holds in the opti- 


mality: 


xij =O and Ina; —Inpy + uj +y% + diz = 0, 


> nig > pyge MI TIAG) > 0, 


or 


xij 20 and Ina; —Inpy + uj +y% + diy = 0, 


—> Liz = pige” Ua tvetnds) = 0. 


The former case holds if p;; = 0; therefore, the optimal solution would be 


written as aj; = pie “7+ +743) > 0 in both cases. 


Note that to obtain the optimal solution of problem (21)—(25), the rela- 
tion obtained for x;; must be considered together with the other optimality 
conditions of the problem, resulting in the following system to be solved: 


Lij = pyge usted) foralliel, ged, 


pe < CY; for all EJ, 

el 

uj(>> Lig — C593) =0 for all j € J, 
ael 

Sy ae for alli € J, 

jEeJ 


n> Ss dizi3 — B) =0, 


jeJ i€l 
dd, titi; SB, 
jeJ ier 
n,u = 0. 


28) 
29) 


30) 
31) 
32) 


33) 


34) 


To transform the nonlinear terms in the system of equations (28)—(34) into 
linear ones, a linearization method is applied. Setting v; =e", s; =e", 
kij = e~"%7 and noting that u,7 > 0, the system of equations (28)—(34) is 


rewritten as follows: 


@yj= pig dj Sikij for alli € I, JE J, 


re < CHU for all jJ€ J; 

iel 

(v; = OD Lig — C593) =0 for all j € J, 
iel 

Stig = W; for alli € J, 


JET 
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(kij -—1I(B-S >So dijxiz)=0 for allie I, je J, (36) 
jeJ i€l 

Laem dita < B, (37) 

jeJ i€l 

vj; <1, kig <1 for allie l,j € J, (38) 


05, 8i,kizg > 0 for alli € I, ged. 


Since u; = O results in v; = 1, for all 7 € J, and visa versa. Therefore, 
constraints (30) and (35) are equivalent. The same statement holds for con- 
straints (32) and (36). Using the approach introduced in [7], the product 
of variables fi; = hijkiz = vjsikiz can be linearized by adding some con- 
straints; see [1]. The final linear mixed-integer system that should be solved 
is equivalent to the following: 


Li = Dig fig for all i € I, ged, (39) 
DE < CY; for all j € J, 
el 
So ay = 0; for alli € J, 
jel 
8, —hiy < Mr;, C593 — >. Big < MA — 1) for ally € J, 
his — fig < Ma, B-SOY dizae <M(1-q) forallieljeJ, 
a 9 
So diay < B, (40) 
jeJ i€l 
fi S hij, hij < 8; for alli € I, j € J, (41) 


Si, haz, fig > 9, 5,0 = {0,1} for alli € J, jE, 


in which M is a sufficiently large number. Problem (17) in Step 1.b can 
be similarly solved while the other problems in Steps 2 and 3 are linear or 
mixed-binary linear that would be solved by existing methods. 


5 Numerical examples 


In this section, the proposed GBD algorithm is implemented in GAM S 23.7 
environment to determine the most probable allocation solution on some 
small and medium-sized networks. During the steps of the algorithm, the 
obtained mixed-binary linear problem (18) is solved by such mixed-integer 
programming solvers as CPLEX. 
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Although by Theorem 1, the algorithm would converge after a finite num- 
ber of iterations, implementing numerical examples revealed that the optimal 
solution would be obtained much sooner than the convergence. Indeed, after 
a number of iterations, the upper bound values might marginally change or 
even remain unchanged, which results in an optimal or near-optimal solu- 
tion. Therefore, we have considered some other stopping criteria, by which 
the algorithm continues until no progress is achieved in the upper or lower 
bounds during a predetermined number of successive iterations, or the total 
number of iterations reaches its maximum value. The used stopping criteria 
are summarized below: 


i. To reach a solution with UB < LB +e, where e€ is a predetermined 
tolerance. 


wt. Proceeding the algorithm with no improvement in the values of UB 
and LB for k successive iterations. 


iti. To reach the maximum number of iterations Mavjter. 


Example 1. Consider the small network shown in Figure 2 wherein 
I =J=VN and the costs of links have been given next to them. The value 
of d;; for each pair (i, 7) is calculated by the shortest distance between i and 
j. The capacities and demands of nodes, (c, w), are inserted in Table 2. To 
estimate the probabilities, the approximate distances with 5 additional cri- 
teria indicated by random integers, representing such geographical and local 
features as proximity to transportation services, security, and surrounding 
facilities, are considered. The compared weights of different nodes j for a 
certain node i with respect to the first attribute, that is, the approximate 
distance, are determined as follows: 


O48 if j ¢ S} = {7 : dij < 25}, 
O85 if j eS? = {j:25 < diy < 35}, 
OF if 7 € 3 = {7:35 < dij < 45}, 
O? if je St = {7:45 < dij < 55}, 
OM if fe SP = {7 +55 < dyh, 


where n? = |S?|, qg = 1,...,5. In addition, the compared weights of nodes 
with respect to the other 5 attributes r = 2,...,6 are given in Table 2. 
We use the multiple attributes decision making procedures to calculate the 
probabilities p;; for different weight vectors W eight, which defines the relative 
importance of 6 criteria; see [1, 29]. 

To see the effect of the available budget as well as emphasizing on dif- 
ferent attributes, the solution of the MIC'pM problem (1)-(6) for p = 3,2 
is estimated for B = 15000 and B = 12000, by applying weight vectors 
Weight, and Weight2; see Table 3. The solutions are estimated via both 
the GBD algorithm with Maxie, = 15, € = 100, and k = 5 as well as the 
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Figure 2: The grid network with 9 nodes 


Table 2: The capacities and demands of nodes along with their weights with 


respect to criteria r = 2,...,6 

[Nodes il 2 3 4 5 6 7 8 9 
ec 300 200 350 320 250 250 350 100 200 

| w 50 60 80 40 90 100 80 60 50 

[| Criteria dl 2 3 4 5 6 7 8 9 
2 0.0541 0.0270 0.1351 0.1622 0.0811 0.1081 0.0541 0.1351 0.2432 
3 0.0244 0.0488 0.1220 0.0976 0.1707 0.1951 0.2195 0.0732 0.0488 
4 0.0455 0.1136 0.1591 0.0227 0.0682 0.1591 0.2045 0.1818 0.0455 
5 0.2162 0.0270 0.0811 0.1351 0.1081 0.0541 0.2432 0.0270 0.1081 
6 0.2286 0.0571 0.1429 0.0286 0.1143 0.0857 0.2571 0.0571 0.0286 


“Dicopt” solver of GAMS 23.7, where the allocation solutions along with 
their objective function values (OBJF’) are inserted in Table 3. 

As it is seen, the optimal values of the objective function provided by the 
GBD algorithm are better than or at least equal to those provided by “Di- 
copt” solver. Note that, the optimal solution might change by changing the 
GAMS solver. However, examining different cases showed that the solutions 
obtained by “Dicopt” are usually better than or at least near to the other 
solvers of GAMS and hence have been reported in Table 3. 

Referring to the estimated solutions, it can be realized if the capacities are 
large enough, the shares of the selected facilities in the MICpM allocation 
solution are proportional to their probabilities. For example, the probabilities 
p(1,9),j =1,...,9 for Weight, are estimated as 


p(1,:) = (0.1617, 0.1072, 0.1090, 0.1423, 0.0983, 0.0895, 0.1471, 0.0642, 0.0808}. 


Among the selected medians 3,6,7, using the GBD algorithm, the median 
7 has the most probability and consequently the most share to provide the 
demand of node 1. Indeed, the higher is the probability p;;, the more would 
be the share of the median v; to provide the demand of node v;. The same 
statement holds for medians 1,3,7 found by the GAMS solver. However, 
decreasing the budget B to 12000 results in a different allocation solution 
to assign as much demand as possible to the closer facilities; compare the 
column of 7 = 6 in cases B = 12000 and B = 15000 for the weight vector 
Weight, in the GBD solution. 

For the weight vector Weight2, wherein the criterion r = 6 is ranked as 
the most important one, the locations of 2 medians are to be selected. As 
it is seen, both the GBD and “Dicopt” solver determine the nodes 3,7 for 
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Table 3: The allocation solution for the MICpM problem 


Facility \ Node T 6 7 OBJIF Time 
B= 15000 p = 3, Weight, = [0.3462 0.1923 0.1154 0.0385 0.1154 0.1923] 
The GBD algorithm 
3 16.739 27.739 36.479 B.3l5 26.951 37.234 12.206 13.266 13.765 2533.350 a 
6 11.952 17.815 29.058 9.391 34.628 40.470 13.927 15.904 18.429 
7 21.308 14.446 14.463 22.294 28.421 22.296 53.867 30.831 17.806 
GAMS solver “Dicopt” 
T 31.095 24.002 15.583 15.006 21.182 16.798 12.873 7.646 5.906 2583.869 3 
3 8.491 27.644 54.607 5.016 36.048 60.544 6.251 12.448 19.226 
9 10.414 8.354 9.810 19.978 32.770 22.659 60.876 39.906 24.869 
B = 12000 p= 3, Weight, = [0.3462 0.1923 0.1154 0.0385 0.1154 0.1923] 
The GBD algorithm 
3 18.245 35.331 46.724 5.194 25.673 37.088 4.160 9.138 12.350 2598.169 a 
6 9.937 17.310 28.391 7.690 43.242 52.847 6.222 14.362 21.675 
7 21.818 7.359 4.885 27.117 21.086 10.065 69.617 36.499 15.975 
GAMS solver “Dicopt” 
T 39.397 34.999 18.405 14.815 14.836 9.744 8.310 4.658 2.983 2655.921 10 
6 2.925 17.626 53.375 4.664 50.319 Trelis 4.155 13.557 26.266 
7 7.678 7.375 20.521 24.845 13.144 67.535 41.785 20.751 
B = 15000 p=2, Weightg = (0.0714, 0.0714, 0.0714, 0.0714, 0.0714, 0.6429] 
The GBD algorithm 
3 20.035 41.620 63.996 7.688 41.899 65.449 6.953 13.305 19.109 2654.990 4 
7 29.965 18.380 16.004 32.312 48.101 34.551 73.047 46.695 30.891 
GAMS solver “Dicopt” 
T 49,154 32.910 5.434 38.559 27.072 6.799 75.123 21.472 3.478 2867.986 5 
3 0.846 27.090 74.566 1.441 62.928 93.201 4.877 38.528 46.522 
B= 12000 P= 2, Weighty =[0.0714, 0.0714, 0.0714, 0.0714, 0.0714, 0.6429] 
The GBD algorithm 
3 21.181 54.045 77.430 3.191 50.387 86.269 0.998 6.688 19.109 2778.247 4 
7 28.819 5.955 2.570 36.809 39.613 13.731 79.002 53.312 30.891 
GAMS solver “Dicopt” 
3 21.181 54.045 77.430 3.191 50.387 86.269 0.998 6.688 19.109 2778.247 4 
id 28.819 5.955 2.570 36.809 39.613 13.731 79.002 53.312 30,891 


B = 12000 as the optimal locations, while for B = 15000, the GBD provides 
a better solution. Referring to Table 2, it is realized that the model finds 
the optimal location of facilities in nodes having more weights with respect 
to the attribute r = 6, while simultaneously the total transportation cost is 
limited to the given budget. In the case of enough capacities, the shares of 
the selected facilities to provide the demands of nodes are again based on 
their corresponding probabilities. In addition, reducing the budget changes 
the allocation solution so that the closer is the median, the larger will be the 
assigned share. 

To analyze the convergence of the GBD algorithm, the variations of 
the upper and lower bounds through successive iterations for weight vec- 
tor Weight, with B = 15000, p = 3 and Weightz with B = 15000, p = 2, 
are illustrated in Figures 3 and 4, respectively. Here, no limit has been con- 
sidered on iterations, that is, the only stopping criterion is satisfying the 
case 7 and the other two ones are neglected. In Figure 3, the algorithm has 
stopped after 76 iterations, where the lower bound becomes greater than the 
upper bound. However, the upper bound values remain almost unchanged 
after 20 iterations. Indeed, the optimal value of the objective function (1) 
has been obtained in the 20th iteration. For the weight vector Weight2 with 
B= 15000, p = 2 in Figure 4, the algorithm converges in the 7the iteration 
while the optimal value of the objective function (1) has been achieved in the 
first iteration. 


Example 2. In the second example, the efficiency of the proposed algorithm 
on larger problems is investigated. First we consider a 100-node network, 
where the network topology, the lengths of links, and the demands and capac- 
ities of nodes are taken from the URL address http : //people.brunel.ac.uk/ 
mastjjb /jeb/orlib/pmedinfo.html. To calculate the probabilities p,;, ap- 
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Figure 3: The variations of the upper and lower bounds for Weight,, B = 
15000, p = 3 


Table 4: Comparison between GBD and GAMS solvers 


B= 800000, p= 10 B=900000, p= 10 B=500000, p= 20 
Solver Objective function Running time | Objective function Running time | Objective function Running time 
Knitro - - - - - - 
Baron - — - - - - 
Bonmin 40639.753 2077 41851.221 1004 40396.948 1006 
Dicopt 41263.215 90 40557.247 16 40340.818 32 
Sbb a a 40247.635 1005 = - 
Lindo global - - - - - - 
Oqulp - - - - - - 
GBD 41929.25 26 41474.553 29 39942.716 26 


proximate distances plus 5 additional criteria, indicated by random integers 
assigned to network nodes, are considered. The problem is solved for p = 10 
and B = 800000, 900000 and p = 20, B = 500000, by applying different 
MINLP solvers of GAMS as well as the GBD algorithm. The stopping 
criteria for the GBD algorithm are the same as those used in Example ?? 
1. The optimal value of the objective function and the elapsed running time 
in each case have been reported in Table 4. The sign (—) indicates that the 
solver has not been able to find any solution for the problem. As it is seen, 
the GBD algorithm generally provides an adequate solution in reasonable 
running time in comparison with the used solvers. 

Next, we use the method proposed by Dolan and More [11] to ana- 
lyze the performance profile of the GBD algorithm as well as three solvers 
Bonmin, Dicopt, and Sbb. To this purpose, we have considered 30 problems 
p © P taken from the URE address http : //people.brunel.ac.uk/ mastjjb 
/jeb/orlib/pmedinfo.html on networks with 100, 250, 300 and 400 nodes, 
and solved them by the three mentioned solvers s = 1,2,3 and the GBD 
algorithm, as well. Let 
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Figure 4: The variations of the upper and lower bounds for Weightz2, B = 
15000, p = 2 


ts,» = the running time required to solve problem p by solver s. 


Referring to [11], the performance on problem p € P by solver s is compared 
with the best performance realized among all the used solvers on this problem. 
Therefore, we use the performance ratio 


Tsp = tewp 
wr Mit fag Pe =e 1 DB) 


A parameter rj > 1s,p for all p,s is chosen, where rs, = rj if and only 
if the solver s does not solve problem p. Then, the probability for the solver 
s that the performance ratio r,,, is less than or equal to a factor Tt € R, is 
demonstrated by 


1 
size{D, Tee < TH. 


ps(T) = |P] 


The function p,(-) is the (cumulative) distribution function for the perfor- 
mance ratio and is termed as the performance profile; see [11]. Figure 5 illus- 
trates the comparison between the performance profiles of the three solvers 
and the GBD algorithm for r = 1,...,200. As it is realized, the GBD algo- 
rithm has much better performance in comparison with the applied GAMS 
solvers. In addition, even for those problems p € P, where the solvers could 
not provide any solution, the GBD algorithm reached a reasonable upper 
bound in less than 200 seconds. 


The most probable allocation solution for the p-median problem 173 


o 
re 


Probability ( Top St) 
oO 
nm 


0.3 


i i i i i n 1 i 
o 20 40 60 80 100 120 140 160 180 200 


Figure 5: Comparison between the performance profiles of the GAMS solvers 
and the GBD algorithm 


6 Summary and conclusions 


In this paper, we proposed a log-based model using the minimum information 
approach, for the capacitated p-median problem to estimate the most likely 
allocation solution, while the total transportation cost should be less than 
or equal to a predetermined budget. Using the generalized benders decom- 
position method, the proposed problem was reduced to easier subproblems 
to which some constraints were added iteratively, and the upper and lower 
bounds have been updated. The convergence of the algorithm was estab- 
lished analytically, while its efficiency in comparison to MINLP solvers was 
investigated practically by applying the performance profile test on some 
numerical examples. 
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